function [stt,etamat,smooth_st,yfor,yferr,logLnc]=kfilterNanSplit_Debug(parvec,funcmod,Y,trainvec,solveopt,addsol)
% ====================================================================
% ====================================================================
%% 1. Model solution
% Matrices of second sample stored in structure second 
[GG1,RR1,C1,eu1,SDX1,ZZ1,~,ssvec,second]=feval(funcmod,parvec,solveopt,addsol);
% A. Non-existence in either sample 
if isequal(eu1,[1;1])==0 || second.euok==0 
    return
end
[Nobs,NY]=size(Y);
%% 2. Demeaning using the split sample 
%
%  First Sample (C1)
if any(C1~=0)==1
    Y(1:addsol.rowBegSecond-1,:)=...
        Y(1:addsol.rowBegSecond-1,:)-repmat((ZZ1*C1)',[addsol.rowBegSecond-1 1]);
end 
% Second sample (second.second.C)
if any(second.C~=0)==1;
    Y(addsol.rowBegSecond:end,:)=...
        Y(addsol.rowBegSecond:end,:)-repmat((second.ZZ*second.C)',[(Nobs-addsol.rowBegSecond+1) 1]);
end 
Y=Y';
%% 3. Forward filter 

%% 3.1. Dimensions and storage 
NZ=size(ZZ1,1); 
NS=size(GG1,1); 
NX=size(SDX1,1); 

vt=zeros(NZ,Nobs);
finvt=zeros(NZ,NZ,Nobs);
kpartg=zeros(NS,NZ,Nobs);
logLnc=zeros(Nobs,1); 
yfor =zeros(NZ,Nobs);

% Matrices with one additional entry (initialization)
% to recover observables 
stt=zeros(NS,Nobs+1);
ptt=zeros(NS,NS,Nobs+1);

W   =eye(NZ);
Zdim      =zeros(Nobs,1); 
mat_obspos=zeros(Nobs,NZ);

%% 3.2 Initialization 
MM=RR1*(SDX1');
stt(:,1)  =zeros(NS,1);
ptt(:,:,1)=lyapunov_symm(GG1,MM*(MM'));


%% 3.3 Start Forward Filter using KF_DK 
% First Part
for ii=1:addsol.rowBegSecond-1;
    
    % Handling of missing observations
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN
    ind =~isnan(ytt);
    rowt =find(~isnan(ytt));
    
    
    ytt=ytt(ind);
    Zdim(ii)=length(ytt);
    Ztt=W((ind==1),:)*ZZ1;
    mat_obspos(ii,1:Zdim(ii))=rowt;
    dimt=( 1:Zdim(ii) );
      
    yfor(dimt,ii) =Ztt*(GG1*stt(:,ii));
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vvv,fin,kgain]...
        =feval(@kf_dk,ytt,Ztt,stt(:,ii),ptt(:,:,ii),GG1,MM);
    
    % Output dimension depends on DIMT 
    vt(dimt,ii)=vvv; 
    finvt(dimt,dimt,ii)=fin; 
    kpartg(:,dimt,ii)=kgain;     
    
    
end
% Change matrix
MM=(second.RR)*(second.SDX');
for ii=addsol.rowBegSecond:Nobs;
    
    % Handling of missing observations
    ytt=Y(:,ii);
    
    % Determine W and position of the NAN
    ind =~isnan(ytt);
    rowt =find(~isnan(ytt));
   
    ytt=ytt(ind);
    Zdim(ii)=length(ytt);
    Ztt=W((ind==1),:)*second.ZZ;
    mat_obspos(ii,1:Zdim(ii))=rowt;
    dimt=( 1:Zdim(ii) );
    
    
    yfor(dimt,ii) = Ztt*(second.GG*stt(:,ii));
    [stt(:,ii+1),ptt(:,:,ii+1),logLnc(ii),vvv,fin,kgain]=...
        feval(@kf_dk,ytt,Ztt,stt(:,ii),ptt(:,:,ii),second.GG,MM);
    
    % Output dimension depends on DIMT 
    vt(dimt,ii)=vvv; 
    finvt(dimt,dimt,ii)=fin; 
    kpartg(:,dimt,ii)=kgain;
    
    
end

%% 4. Likelihood with Integration Constant  
consInt=-0.5*sum(Zdim)*log(2*pi); 
logLnc=logLnc(trainvec(1):trainvec(2)); 
logL=consInt+sum(logLnc);

%% 5. Truncate filters and obtain initial observations 
yferr=vt';
yfor =yfor'; 
azero=stt(:,1);Pzero=ptt(:,:,1);
stt=stt(:,2:end); ptt=ptt(:,:,2:end); 

%% 6. Disturbance smoother with TV matrices 
% Obtain the Innovations using a disturbance smoother 
etamat=zeros(NX,Nobs); 
smooth_st=zeros(NS,Nobs); 
rmat=zeros(NS,Nobs); 

Gf=second.GG;
Zf=second.ZZ;
Rf=second.RR;
Qf=(second.SDX')*second.SDX;

%% 6.1 Initialize RSTAR & start at t=Nobs 
rstar=zeros(NS,1); 
[rstar,etamat(:,end)]=smoothdis(rstar,Qf,(Rf'),(Zf'),finvt(dimt,dimt,end),zeros(NS),vt(dimt,end));
smooth_st(:,ii)=stt(:,end)+ptt(:,:,end)*rstar; 
rmat(:,end)=rstar; 

%% 6.2 Begin Backward recursion 
Ztr=Zf'; 
for ii=(Nobs-1):-1:1;
    
    %% Volatility must always be at time t 
    if ii==addsol.rowBegSecond-1;
        Ztr=ZZ1';
        Rf=RR1;
        Qf=(SDX1'*SDX1);
    end
    
    % Varying dimension in the backward recurions
    dimt=( 1:Zdim(ii) );
    % [Ztt]'= [Wtt*Z]' = = Z'*Wtt'
    Ztt=Ztr*( W( mat_obspos(ii,1:Zdim(ii)),:)');
    
    
    [rstar,etamat(:,ii)]=smoothdis(rstar,Qf,(Rf'),...
        Ztt,finvt(dimt,dimt,ii),((Gf-Gf*kpartg(:,dimt,ii)*Ztt')'),vt(dimt,ii));
    smooth_st(:,ii)=stt(:,ii)+ptt(:,:,ii)*rstar;
    rmat(:,ii)=rstar; 
    % Recall addsol.rowBegSecond -1 is the end of the First Sample
    % This break must occur exactly when first sample ends
    if ii==addsol.rowBegSecond-1;
        Gf=GG1;
    end
end
etamat=(etamat)';
smooth_st=(smooth_st)'; 
stt  =stt'; 


%% 7. Check Smoother 
% Check that Smooth States are identical if using disturbance smoother (above) vs. state smoother (below) 
% and if can also recover the observables
outcount.aold=azero;
outcount.Pold=Pzero;
tol=1e-5; 
% State at time zero give by GG1*azero+Pzero*rstar;
%% $$ s(1)=GG*s(0) + n(1) $$ where $$ n(1)=P(0)*r(1) $$ 
%% $$ s(0|T)=s(0|0) + P(0|0)G(1)'r(1) $$ 
%% $$ s(1|T)=GG*s(0|T) + R(1)n(1|T) $$
[yCheck,smoothCheck]=kfilterRegSplitSimulation(azero+Pzero*GG1'*rstar,etamat'); 
disp('Max Discrepancy Smooth vs. Actual Data'); 
maxdifY=comparemat(yCheck,Y'); 
disp('Max Discrepancy State and Innovation Smoother'); 
maxdifS=comparemat(smoothCheck,smooth_st); 
if maxdifY > tol || maxdifS > tol 
    error('Smoother discrepancy exceeds tolerance') 
end 

%% 8. Initial Variance and State per shock (used below for the counterfactual decompositions)
% vInitialPershock: Initial Variance, decomposed per shock 
% sOnePerShock, State at time 1 decomposed, decomposed per shock
vInitialPerShock=zeros(NS,NS,NX);
for ii=1:NX
    vInitialPerShock(:,:,ii)=lyapunov_symm(GG1,RR1(:,ii)*SDX1(ii,ii)*SDX1(ii,ii)*(RR1(:,ii)'));
end
sOnePerShock=zeros(NS,NX);
for ii=1:NX
    sOnePerShock(:,ii)=vInitialPerShock(:,:,ii)*rstar;
end
maxdifSzero=comparemat(sum(sOnePerShock,2),smooth_st(1,:)' ); 
if maxdifSzero > tol; error('Cannot decompose sonePerShock'); end 

%% 9. Counterfactual Decomposition 
% Generate states by feeding each shock at a time, ensure that it recovers
% the original state 
disp('Begin Counterfactual')
countStates=zeros(Nobs,NS,NX);
countObs   =zeros(Nobs,NY,NX);
for ii=1:NX 
    etaTemp=zeros(Nobs,NX); 
    etaTemp(:,ii)=etamat(:,ii); 
    [countObs(:,:,ii),countStates(:,:,ii)]=kfilterRegSplitSimulation(sOnePerShock(:,ii),etaTemp');
end 
disp('Maximum Discrepancy Counterfactual States and Smooth States') 
maxdifCount=comparemat(sum(countStates,3),smooth_st); 
if maxdifCount > tol; error('Counterfactual States do not recover smooth states'); end 
    
   
%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [NS 1] Initial State 
% innovMat: [Nx Nobx] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat) 
        % Note: Rememeber that feeding in S(1) not S(0) where 
        % S(1)=GG*S(0)+P(0)*r(1)
        if ~isequal(size(innovMat),[NX Nobs]) 
            error('Input innovMat must be [NX Nobs]') 
        end 
        sSim=zeros(NS,Nobs);
        ySim=zeros(size(Y));
        sSim(:,1)=GG1*sInitial+RR1*innovMat(:,1); 
        ySim(:,1)=ZZ1*sSim(:,1);
        for kk=2:addsol.rowBegSecond-1;
            sSim(:,kk)=GG1*sSim(:,kk-1)+RR1*innovMat(:,kk);
            ySim(:,kk)=ZZ1*sSim(:,kk);
        end
        for kk=addsol.rowBegSecond:Nobs;
            sSim(:,kk)=second.GG*sSim(:,kk-1)+second.RR*(innovMat(:,kk));
            ySim(:,kk)=second.ZZ*sSim(:,kk);
        end       
        ySim=ySim'; 
        sSim=sSim'; 
    end
end 